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ABSTRACT 

Subsonic  and  transonic  steady  and  unsteady  flowfields  over  airfoils  are  inves- 
tigated with  the  numerical  solution  of  the  governing  equations.  This  study  aims  to 
enhance  the  performance  of  rotary  wing  and  fixed  wing  aircraft  by  better  understand- 
ing and  by  taking  advantage  of  unsteady  phenomena  such  as  dynamic  lift.  In  the  past 
few  years  many  advances  have  been  made  in  algorithm  development  for  the  numerical 
solution  of  the  Euler  and  the  Navier  Stokes  equations.  In  this  study,  these  new  zonal 
techniques  axe  applied.  A  zonal  approach  is  more  computationally  efficient  in  solving 
the  governing  equations  than  previous  approaches,  and  has  certain  advantages  over 
the  standard  single  moving  grid  approach.  The  zonal  grid  consists  of  two  grids,  one 
being  the  inner  grid  which  is  fixed  to  the  airfoil,  and  the  other  being  the  outer  grid 
which  extends  to  the  far  field  or  to  a  specified  outer  boundary.  The  inner  grid  is 
allowed  to  rotate  with  the  body,  while  the  outer  grid  remains  fixed.  The  thin-layer 
Navier-Stokes  equations  are  solved  for  the  inner  grid,  while  the  Euler  equations  are 
solved  for  the  outer  grid.  Communication  between  the  two  grids  is  accomplished  by 
interpolating  the  flow  quantities  at  the  zonal  interface.  Solutions  axe  obtained  for 
flows  at  fixed  angles  of  incidence,  and  for  unsteady  flows  over  pitching  and  oscillating 
airfoils.  The  computed  results  are  in  good  agreement  with  available  experimental 
data. 
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I.  INTRODUCTION 

A.   BACKGROUND 

Investigation  of  steady  and  unsteady  flowfields  over  airfoils  is  an  active  area 
of  numerical  and  experimental  research.  Unsteady  pitch-up  motion  of  airfoils  alters 
significantly  the  aerodynamic  characteristics  of  lifting  surfaces.  Pitch-up  motion  of 
an  airfoil  produces  lift  augmentation  and  delay  of  stall  to  higher  angles  of  incidence 
compared  to  airfoils  held  at  a  steady  incidence.  Understanding  the  mechanisms  that 
cause  the  dynamic  lift  development  is  a  subject  of  interest  in  both  theoretical  and 
applied  research. 

A  comprehensive  explanation  of  the  dynamic  lift  phenomenon  is  given  by  Tyler 
and  Leishman  [Ref.  lj.  To  briefly  summarize,  unsteady  airfoil  motion  tends  to 
maintain  high  lift  to  higher  angles-of-attack,  because  the  unsteady  flow  causes  a 
time  delay  in  the  build-up  of  the  lift  force  and  the  adverse  pressure  gradient.  The 
unsteady  motion  also  gives  the  airfoil  a  virtual  camber  which  decreases  the  leading 
edge  pressure  and  pressure  gradients.  Depending  on  the  flow  conditions,  as  the  airfoil 
begins  to  stall,  a  vortex  forms  at  the  leading  edge  which  grows  with  time  which  is 
eventually  shed  downstream  in  the  wake.  This  vortex  shedding  phenomenon  alters  the 
chordwise  pressure  distribution  on  the  upper  surface  of  the  airfoil  resulting  in  higher 
maximum  lift  coefficients.  Sometimes,  however,  only  a  recirculating  flow  is  observed 
at  the  airfoil  trailing  edge.  As  the  pitch-up  motion  progresses  the  recirculating  region 
extends  upstream  towards  the  leading  edge  and  eventually  a  vortex  is  formed  at  the 
trailing  edge. 


Rotary  wing  and  fixed  wing  aircraft  designers  can  enhance  the  performance  of 
the  aircraft  by  taking  advantage  of  the  dynamic  lift  concept.  However  ,  the  resulting 
undesirable  pitching  moment  variations  have  deprived  helicopters  and  airplanes  of 
the  benefits  of  dynamic  lift  [Ref.  2].  Carr  [Ref.  3]  gives  a  comprehensive  review  of 
the  progress  that  has  been  made  in  the  study  of  the  dynamic  stall  phenomenon. 

Experimental  work  on  steady  transonic  flow  on  airfoils  was  conducted  by  McDe- 
vitt  and  Okuno  [Ref.  4].  These  data  cover  a  wide  range  of  flow  conditions  and  can 
be  used  for  steady  code  validation.  One  of  the  pioneering  experimental  works  on  un- 
steady airfoil  flows  was  conducted  by  McCroskey  et  al  [Ref.  5].  In  their  experiments, 
unsteady  dynamic  effects  for  two-dimensional  airfoil  flows  were  studied  for  the  first 
time  in  conjunction  with  flows  over  helicopter  blades.  They  showed  the  effects  of 
unsteady  lift  and  pitching  moment  on  the  retreating  rotor  blades.  Their  research 
also  showed  that  unsteady  motion  parameters  such  as  reduced  frequency,  appeared 
to  be  more  important  than  airfoil  shape  in  determining  the  dynamic  stall  airloads. 
Benchmark  experimental  data  of  oscillating  and  pitching  airfoils  was  also  collected  by 
Landon  [Ref.  6].  The  experimental  data  of  references  [5]  and  [6]  gives  enough  quali- 
tative information  for  unsteady  code  validation.  Recent  experiments  with  oscillating 
airfoils,  performed  by  Chandrasekhara  and  Brydges  [Ref.  7],  have  shown  definitively 
that  at  subsonic  Mach  numbers  the  unsteady  flow  in  the  vicinity  of  the  leading  edge 
can  reach  supersonic  speeds  and  generate  shocks. 

Along  with  all  the  experimental  studies  being  conducted,  a  great  effort  is  also 
underway  to  compute  unsteady  viscous  flows.  Developments  of  numerical  methods 
for  the  Navier-Stokes  equations  [Refs.  8,  9,  10]  during  the  past  few  years  provide 
new  tools  for  the  investigation  and  prediction  of  airfoil  flows.  In  this  investigation 
the  compressible  thin  layer  Navier-Stokes  equations  are  solved  using  a  zonal  grid 
approach.    The  objective  is  to  develop  a  computationally  efficient  method  to  study 


two-dimensional  unsteady  flows.  This  will  be  accomplished  by  developing  a  solution 
procedure  for  two  grids,  an  inner  viscous  grid  around  the  solid  body,  and  an  outer 
grid  which  is  coarser  representing  the  outer  flow  field.  The  inner  grid  is  free  to  rotate 
to  any  angle  of  attack,  and  the  outer  grid  remains  stationary. 

The  idea  of  moving  meshes  and  dynamic  meshes  is  not  new.  The  dynamic 
and  adaptive  grid  solution  method  refers  to  computational  grids  which  are  coupled 
to  the  physical  problem  that  is  being  solved. An  adaptive  solution  procedure  with 
grid  points  that  continually  move  during  the  solution  process  in  order  to  resolve  the 
developing  gradients  has  been  shown  [Ref.  11].  A  dynamic  type  of  mesh  has  been 
applied  by  Rumsey  and  Anderson  [Ref.  12]  to  simulate  aileron  buzz  using  the  thin 
layer  Navier-Stokes  equations. 

The  idea  of  overlapped  or  patched  grid  schemes  was  used  by  Rubbert  and  Lee 
[Ref.  13]  with  the  limitation  that  the  grid  lines  at  the  boundaries  were  continu- 
ous. Benek  et  al  [Ref.  14]  developed  the  "CHIMERA"  approach  for  two-dimensional 
embedded  and  overlapping  grids.  They  demonstrated  that  the  grid  lines  at  the  over- 
lapped boundaries  did  not  have  to  be  continuous,  and  that  the  flow  quantities  could 
be  successfully  interpolated.  The  "CHIMERA"  approach  was  later  extended  to  three 
dimensions  by  Benek  et  al  [Ref.  15]  .  Rai  [Ref.  16]  developed  a  technique  for  in- 
dependent zonal  grids  where  the  flow  variables  were  interpolated  across  the  zonal 
interface.  Chesshire  and  Henshaw  [Ref.  17]  developed  a  methodology  for  solving  the 
steady  compressible  Navier-Stokes  equations  using  multiple  overlapped  grids.  Their 
methodology  was  demonstrated  by  solving  the  governing  equations  on  a  composite 
grid  that  was  comprised  of  an  airfoil  grid,  a  leading-edge  flap  grid,  a  trailing-edge 
flap  grid  and  a  grid  for  the  outer  flowfield.  Chyu  and  Davis  [Ref.  IS]  investigated 
unsteady  transonic  flow  using  a  moving  airfoil  grid.  They  developed  stationary  com- 
putational grids  around  the  airfoil  at  its  lowest  and  highest  angles  of  attack,  and  then 


interpolated  a  new  grid  as  the  airfoil  oscillated  to  intermediate  angles-of-attack.  Reu 
and  Ying  [Ref.  19]  developed  a  composite  grid  approach  to  study  the  flow  about 
pitching  airfoils  in  a  wind  tunnel.  Their  approach  consisted  of  a  structured  inner  grid 
and  an  unstructured  grid  for  the  outer  flowfield. 

Unsteady  problems  have  also  been  solved  by  oscillating  the  incoming  flow,  as 
opposed  to  moving  the  grid.  However,  this  method  does  not  work  if  a  solution  is 
sought  for  two  or  more  objects,  such  as  a  wing  tail  combination  or  a  canard  wing  in 
relative  motion  to  each  other. 

B.      PURPOSE 

In  this  work,  unsteady  compressible  flows  are  investigated  using  a  numerical 
technique  which  is  applied  to  zonal  meshes.  The  governing  equations  are  solved 
on  multiple  computational  grids,  where  one  of  the  grids  is  free  to  move  in  unison 
with  the  solid  boundary  and  the  other  grid  is  fixed.  The  meshes  overlap  at  the 
zonal  interface.  The  scheme  that  is  developed  is  similar  to  the  "CHIMERA"  scheme 
mentioned  previously,  but  we  impose  the  restriction  of  a  circular  shape  on  the  zonal 
interface. 

The  new  approach  developed  in  this  study  is  more  computationally  efficient 
compared  to  previous  schemes  used  to  study  unsteady  flows.  The  zonal  grid  approach 
avoids  the  need  to  regenerate  or  interpolate  entire  grids  at  every  angle  of  attack.  The 
zonal  interface  is  in  a  known  position;  therefore  the  entire  flowfield  does  not  need  to 
be  searched.  In  addition,  since  the  outer  grid  remains  stationary,  the  metric  terms  do 
not  need  to  be  recomputed  at  each  time  step.  Another  advantage  is  that  the  zonal 
interface  is  circular;  therefore  the  flow  variables  are  interpolated  in  the  circumferential 
direction  only.  The  zonal  grid  approach  also  allows  for  the  application  of  different 
solution  methods  to  the  inner  and  outer  grids.  For  example,  viscous  solution  methods 


near  the  solid  body  and  an  inviscid  method  on  the  outer  mesh  can  be  implemented. 

The  primary  objective  of  this  investigation  is  to  develop  and  test  the  zonal  grid 
methodology.  The  space  discretization  is  based  on  Osher's  [Ref.  20]  upwind  method. 
An  implicit  scheme  is  used  for  time  integration.  An  advantage  of  upwind  schemes  is 
that  they  are  naturally  dissipative  and  no  explicit  artificial  dissipation  is  required. 

In  order  to  meet  our  goal,  a  procedure  for  generating  zonal  grids  was  devel- 
oped, along  with  the  zonal  flow  solver.  The  effect  of  grid  resolution  on  the  accuracy 
of  solutions  was  studied  first  for  steady  flows  and  the  solutions  were  compared  to 
experimental  data  collected  by  Harris  [Ref.  21].  The  unsteady  results  were  validated 
using  the  experiments  by  Landon.  The  dependence  of  the  solution  on  the  location  of 
the  overlapped  zonal  region  was  examined.  Viscous  steady  state  solutions  were  com- 
puted and  compared  to  experimental  data.  Finally,  unsteady  flows  were  investigated 
by  computing  solutions  for  airfoils  in  ramp  and  oscillatory  motion.  The  ramp  motion 
was  started  at  zero  angle  of  incidence,  and  ramped  up  to  30  degrees.  The  computed 
pressure  coefficients  were  verified  by  comparing  with  available  experimental  pressure 
distributions.  The  boundary  layers  computed  for  this  case  were  also  compared  to  an 
interactive  boundary  layer  code  [Ref.  22].  The  oscillatory  test  case  was  verified  by 
comparing  the  computed  pressure  coefficients  with  experimental  data.  This  approach 
and  the  computed  results  are  described  in  the  following  chapters. 


II.  GOVERNING  EQUATIONS 

In  order  to  compute  compressible  viscous  fluid  flow  around  a  body,  the  conti- 
nuity, momentum  and  energy  equations  must  be  solved  simultaneously.  The  vector 
and  the  conservation-law  form  of  the  compressible  Reynolds-averaged  Navier-Stokes 
equation  is  presented.  A  detailed  derivation  can  be  found  in  [Ref.  23] 

A.  CONTINUITY  EQUATION 

The  continuity  equation  expresses  the  conservation  of  mass  law  applied  to  a 
fluid  passing  through  a  control  volume  fixed  in  space 

^  +  (V-pV)=0  (2.1) 

here  p  is  the  fluid  density  and  V  is  the  fluid  velocity.  Equation  2.1  states  that  the 
net  mass  flux  through  a  control  volume  bounding  surface  must  be  equal  to  the  time 
rate  of  change  of  the  mass  inside  the  control  volume.  In  a  two-dimensional  Cartesian 
coordinate  system  this  equation  reads 

%  +  |-(H  +  ^(pw)  =  0  (2.2) 

at      ox  oz 

where  u  and  w  are  velocity  components  along  the  x  and  z  directions,  respectively. 

B.  MOMENTUM  EQUATION 

The  momentum  equation  expresses  Newton's  second  law  as  applied  to  a  fluid 
element  passing  through  a  control  volume  fixed  in  space.  The  momentum  equation 
is: 

-(pV)  +  V  •  pVY  =  pt  +  V  •  Ily  (2.3) 


The  first  term  in  equation  2.3  represents  the  time  rate  of  change  of  momentum 
per  unit  volume  in  the  control  volume.  The  second  term  represents  the  moment  flux 
through  the  bounding  surface  of  the  control  volume.  /  is  the  body  force  per  unit 
volume  and  Uv  is  the  stress  tensor  given  by 


n,j  =  ~p$ij  +  H 


du{      2      duk 

—  -Oi 


dxj       3  tJdxk 
where  i,j,k  —  1,2,3  and  StJ  is  the  Kronecker  delta. 


(2.4) 


By  substituting  equation  2.4  into  equation  2.3  and  expanding  equation  2.3  foi 
a  two-dimensional  Cartesian  coordinate  svstem  obtain 
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nDw    —    nf  ^2  _L  JL    \ln  f^^w  9u\]     1     _d_   \      (  dm     ,     du\ 

P~dT  ~  PJ=  ~  tz  +  Tz  [3P  [-JI  ~  Tx)\+  Tx  [P  {-dx-  +  Tz). 
Equations    2.5  are  known  as  the  Navier-Stokes  equations  for  two-dimensional  flow. 

C.      ENERGY  EQUATION 

The  energy  equation  is  derived  by  applying  the  first  law  of  thermodynamics  (  rate 
of  change  of  energy  =  net  heat  flux  into  particle  -f  rate  of  work  done  on  particle.   ) 

de      dQ         .  .  ,     .  d  .  _  . 

-77  -  -z-  -  p[fxU  +  JzW)     +     tt- (eu  +  pu  -  utxx  -  wrrz  +  Qx) 
at       at  ox 

+     —(ew  +  pw  —  wtZz  —  uTj.-  +  Q,)  =  Q 

where  e  is  the  total  energy  per  unit  volume  and  Q  is  the  heat  addition  per  unit 
volume. 

The  above  equations  can  be  rewritten  in  non-dimensionalized  vector  form  as 
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dt        dx        dz        Re  \  dx 
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where  Uoo  is  the  free  stream  speed  and  L  is  a  reference  length  The  pressure  is  related 
to  u,  w,  p,  and  e  by 

p  =  (7  -  l)[e  -  -p(u2  +  w2)] 

In  the  above  equations,  the  ratio  of  the  specific  heats,  7,  is  set  equal  to  1.4  and 
a2  —  ~ip/p  is  the  local  speed  of  sound. 

The  density  is  non-dimensionalized  by  the  free  stream  density,  /?oo,  the  velocities 
are  non  dimensionalized  by  the  free  stream  speed  of  sound,  and  the  total  energy  is 
non  dimensionalized  by  Pooalo- 
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D.      TURBULENCE  MODEL 

The  Navier-Stokes  equations  can  completely  model  fluid  flow,  however,  in  or- 
der to  resolve  turbulent  scales  at  high  Reynolds  numbers  and  realistic  geometries. 
very  high  grid  densities  are  required.  Therefore,  present  state-of-the-art  algorithms 
and  computer  technology  allows  direct  simulation  of  turbulent  flows  for  only  simple 
geometries  and  low  Reynolds  number  flows  which  are  of  limited  practical  interest. 
In  order  to  enable  computation  of  turbulent  flows  for  configurations,  of  practical  in- 
terest, turbulence  modeling  is  used.  Turbulence  models  are  implemented  with  the 
time-averaged  forms  of  the  Navier-Stokes  equations. 

Two  widely  used  averaging  procedures  are  the  standard  time  averaging  proce- 
dure for  incompressible  flow,  and  the  mass-averaged  approach  for  compressible  flows 
[Ref.  24].  The  time  averaging  procedure  destroys  high  frequency  information  of  the 
turbulence,  but  the  unsteady  mean  flow  information  is  preserved. 

For  the  incompressible  case,  the  randomly  changing  flow  variables  are  replaced 
with  their  averages  plus  their  fluctuations.  For  example,  in  the  Cartesian  coordinate 
system,  the  u  velocity  component  is  represented  as  u  =  u  -f-  u',  where  u  is  the  mean 
velocity  and  u'  is  the  fluctuation  about  the  mean.  The  governing  equations  are  time 
averaged,  and  the  average  of  the  fluctuation  terms  is  set  equal  to  zero. 

In  the  compressible  flow  case,  the  mass-weighted  variable  of  the  Favre  averaging 
approach  is  used.  In  this  case,  u  is  represented  as  u  =  u  +  u"  where  u  is 

pu 

P 

Here  the  time  average  of  the  doubly  primed  fluctuating  quantities  is  not  equal  to  zero. 
After  the  substitution  is  carried  out  for  all  of  the  fluctuating  flow  variables,  the  entire 
equation  is  time  averaged.  Next,  all  the  time  averaged  terms  that  are  doubly  primed 


and  multiplied  by  density  are  defined  to  be  zero.  For  example 

pu"  =  0 

The  equations  of  mean  motion  resulting  from  the  time  averaging  procedure  have 
more  unknowns  than  equations.  This  constitutes  the  closure  problem  of  turbulence. 
In  order  to  close  these  equations  a  turbulence  model  is  used.  Several  models  have 
been  proposed,  in  this  study. 

The  Baldwin-Lomax  (B  —  L)  turbulence  model  [Ref.  25]  was  used.  This  model 
is  a  two-layer  eddy  viscosity  model  which  simulates  the  effect  of  turbulence  in  terms 
of  the  eddy  viscosity  coefficient  \it.  The  term  \i  in  the  stress  terms  is  replaced  by 
(i  +  fit,  and  the  f.i/Pr  in  the  heat  flux  terms  is  replaced  with  /i/  Pr  +  [it/  Prt.  The  B  —  L 
turbulence  model  is  similar  to  the  Cebeci-Smith  turbulence  model  [Ref.  26],  but  it 
bypasses  the  need  for  finding  the  edge  of  the  boundary  layer  by  using  vorticity  instead 
of  the  boundary  layer  thickness.  This  model  is  adequate  for  flows  which  have  mild 
pressure  gradients,  but  it  is  not  very  suitable  for  highly  separated  flows.  A  complete 
description  of  the  model  is  given  in  reference  [25].  The  basic  equations  of  the  model 
follow. 

In  the  inner  layer,  the  eddy  viscosity  is  assumed  to  be  proportional  to  the  mixing 
length  and  vorticity,  and  in  the  outer  layer  it  uses  an  exponentially  decaying  formula. 
The  inner  eddy  viscosity  is  computed  up  to  the  point  where  it  is  equal  to  the  outer 
eddy  viscosity  as  shown  below. 

{l/'t/inner       V   2z   ^/crossover  I  O  "  \ 

(f-^t  Jouter      2/crossover  ^  V 

where  y  is  the  normal  distance  from  the  wall  and  ^crossover  taken  at  its  minimum  value 
where  it  equals  y.  The  inner  eddy  viscosity  is  given  by: 

(^t)inner  =  P^M 
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wnere 


/  =   Klj 


1  —  exp(- 


.4+ 


u.' 


N 


dw      du 
dx       dz 


+  PwUrV 

y   = 


,4+  is  an  experimentally  determined  damping  constant,  k  is  the  Von  Karman  constant 
The  outer  eddy  viscosity  is  given  by 


/'outer  =  KCcppF(y)\VAKEF{y)KLEB- 


F(y)KLEB  is  the  Klebanoff  intermittency  factor  given  by 


F(y) 


KLEB 


1  +  5^(CklebI 


-I  -1 


«  and  Ccp  are  constants.  For  boundary  layers 


"  \y  J  wake  —  J/max-*  max- 


For  wakes  and  separated  boundary  layers 


i?(y)wake  =  Cwkyt 


lDlF 


Fn 


The  quantity  ym&x  is  the  value  of  y  determined  for  the  maximum  value  of  Fmax  and 
Fmax  is  determined  by 


F(y)  =  y\ 


1  —  exp    -4+ 


E.      TRANSFORMATION  TO  GENERALIZED  COORDINATES 

In  order  to  use  an  unweighted  differencing  scheme  that  facilitate  the  numerical 
implementation  on  body  fitted  coordinate  systems  suitable  for  complex  configura- 
tions, the  equations  are  transformed  to  a  generalized  coordinate  system  using  the 
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following  transformations 


£  =  £(x,z);      C  =  ((x,z). 


[2.8) 


The  above  equations  transform  the  governing  equations  from  the  physical  domain  to 
a  body  fitted  coordinate  system.  The  transformation  is  carried  out  by  using  the  chain 
rule  of  partial  differentiation 


d_     f  d_     .  d_      d_        d_     .  d_ 

dx  ~  sx  cX  +  C*  dQ ;      dz  ~  iz  d£  +  G  d( 


(2.9) 


For  example,  the  continuity  equation  would  be  transformed  in  the  following  manner. 


dpu  dpu       ,  dpv       dpv  dpu       „  dpv 

~dx~  =  'x~dl  +  Cx^cr;    ~d7  =  ^~dl  +  ^~ec 


[2.10) 


(r,  £z  Xx  and  (\  are  known  as  the  metrics.  The  metric  terms  can  determined  in  the 
following  manner.  First  write  the  differential  expressions  for  £  and  ( 


d£  =  ^xdx  +  £zdz]      d(  =  (xdx  +  (zdz 


2.11) 


In  matrix  form  the  expressions  become 


d( 


6    6 


Next  we  write  the  differential  expression  for  x  and  z 


'2.12: 


dx  =  x^d^  +  x^c/£ 


In  matrix  form  thev  become 


dx 

dz 


d( 
dC 


By  comparing  equations  2.12  and  2.14  we  can  write 

-l 


XZ     x< 


6    6 


2.13: 


(2.14; 


,2.15) 
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and  solving  for  the  metrics  we  get 


where 


J 


Sj         —   J  -c, 
Cx       =  ~Jzt 

G     =  Jzt 

1 


/2.16) 


xz    xc 

The  mapping  from  (x,  z)  to  (£,  (,')  is  one  to  one  if  the  Jacobian,  J,  is  non  singular. 


(2.17) 


F.      THIN  LAYER  NAVIER  STOKES  EQUATIONS 

In  order  to  accurately  calculate  the  normal  gradients  in  the  boundary  layer,  it 
is  necessary  to  make  the  normal  grid  spacing  very  fine  close  to  the  solid  surfaces,  n 
the  streamwise  direction  the  flow  gradients  are  not  large  and  no  fine  grid  spacing  is 
required.  As  a  result  the  grid  cells  near  the  body  have  a  very  high  aspect  ratio.  With 
this  type  of  grid,  existing  gradients  in  the  streamwise  direction  are  not  fully  resolved. 
In  order  to  facilitate  the  numerical  implementation  of  the  thin  layer  approximation 
is  employed  by  retaining  only  the  viscous  derivatives  normal  to  the  body.  The  thin 
layer  Navier-Stokes  equations  transformed  into  a  generalized  curvilinear  system  is 
vector  form  are 


where 


dQ      OF      dG 


Q  = 


1 

Re 

dS 

P 

1 

pu 

J 

pw 

e 

F  = 


PU 

pull  +  £xp 

pwU  +  (xp 

[e  +  p)U-£tp 
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G  - 


pW 

puW  +  (xp 
pwW  +  (zp 

_  {e+p)W  -Qp  . 

The  viscous  flux  term  is  transformed  as 


1 

7 


s  = 


1 
1 


w 


here 


0 

j.imiuc  +  {fi/3)m2(x 

firriiw^  +  (^/3)m2Cx 

prrtim^  +  (^/3)m2m4 


mi  =  c  +  c2 


m2  =  Cx^c  +  (z^C 


3  2  3£  2  "l"    Pr(-y-l)    ac 


(2.18; 


m4  =  ("xw  +  C2u;. 
U  and  VK  are  the  contravariant  velocity  components.  These  velocity  components 

are  normal  and  parallel  to  the  constant  £  and  (  surfaces,  respectively. 

U  =  6  +  6u  +  k"> 

W  =  C«  +  uw  +  C*w 
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III.  SOLUTION  METHODS 

A.      ZONAL  GRID  GENERATION 

The  zonal  grids  were  generated  using  a  software  package  called  GRAPE2D.  de- 
veloped at  NASA  Ames  Research  Center  by  Reese  L.  Sorenson  [Ref.  27].  GRAPE2D 
generates  field  grids  by  solving  the  following  set  of  elliptic  Poisson  equations: 

U~U  =  n  (3-1) 

Vxx   -  Vzz    =   Q 

GRAPE2D  was  used  to  generate  the  inner  and  outer  grid  separately.  In  gener- 
ating the  inner  grid,  the  number  of  circumferential  and  radial  grid  points,  the  spacing 
at  the  body  surface,  and  the  radius  and  shape  of  the  outer  boundary  were  specified. 
Once  the  inner  grid  was  generated,  the  outer  boundary  of  the  inner  grid  was  used  as 
the  inner  boundary  for  the  outer  grid. 

Due  to  the  placement  of  the  reference  axis  at  the  quarter  chord  of  the  airfoil, 
the  grid  spacing  in  front  of  the  airfoil  is  larger  than  aft  of  the  airfoil.  This  causes  a 
problem  in  specifying  an  initial  grid  spacing  for  the  outer  grid.  In  order  to  obtain 
the  smoothest  overall  transition  from  the  inner  grid  to  the  outer  grid,  the  starting 
spacing  for  the  outer  grid  was  obtained  by  averaging  the  minimum  and  maximum 
spacing  at  the  outer  boundary  of  the  inner  grid.  The  outer  boundary  of  the  outer 
grid  was  prescribed  as  a  rectangular  region,  usually  about  six  chord  lengths  away 
from  the  body. 

The  perfectly  circular  zonal  boundary  was  accomplished  by  reading  the  over- 
lapped regions  x  and  z  coordinates  and  finding  an  average  radius  which  was  measured 
from  the  center  of  the  grid,  regardless  of  the  location  of  the  inner  body.  Next,  an 
arbitrary  grid  point  was  chosen  to  use  as  a  reference  point.    Then  an  angle  theta, 
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0,  was  calculated  which  represented  the  angle  between  the  reference  grid  point  and 
the  grid  point  under  consideration.  Then  having  the  average  radius  and  the  angle  0 
for  the  grid  points  in  the  overlapped  region,  a  new  set  of  x  and  y  coordinates  was 
calculated  using  the  following  equations: 


x  =  r  cos  6  ,  y  —  r  sin  0 


(3.2 


The  overlapping  of  the  two  grids  was  performed  by  adding  the  next  to  last  layer 
of  grid  points  from  the  inner  grid  to  the  outer  grid,  and  adding  the  second  layer  of 
grid  points  from  the  outer  grid  to  the  outer  boundary  of  the  inner  grid.  This  resulted 
in  an  overlapped  region  of  three  grid  points  or  two  grid  cells  in  the  radial  direction. 
The  grid  was  also  overlapped  by  a  grid  point  in  the  circumferential  direction. 

B.     NUMERICAL  SCHEME 

The  numerical  integration  is  performed  using  an  upwind-biased,  factorized. 
iterative,  implicit  numerical  scheme  [Ref.  20]  given  by 

[/  +  MV^  +  A{A-fc)]Px 

-[(Qlk-Qlk)  +  hdF[+l/2,k-^-i/2,k) 
+Ac(G?,*+i/2  -  <5?.*-i/J  "  R^hc(Sl+1/2  -  Sl_l/2)} 

In  equation  3.3,   h^   =   Ar/A£,  etc.,  A±   =   (dF/dQ),  etc.,  are  the  flux  Jacobian 

matrices,  and  A,  V,  and  8  are  the  forward,  backward  and  central  difference  operators, 

respectively.  The  quantities  F,+i/2,;-,  Giik+i/2,  and  Sitk+i/2  are  numerical  fluxes.  The 

inviscid  fluxes  F  and  G  are  evaluated  using  Osher's  upwinding  scheme  [Ref.  19].  The 

numerical  fluxes  for  a  third-order  accurate  scheme  are  given  by 


Fi+i/2,k  =  Fi+i/2,k  +  |  [AFtt1/2)k  +  2^i++i/2./fc 
I   ^;3/2l*  +  2AFr+1/2, J  =  F{Qitk,Qi+i,k)  + 
|  [AF+(Qi+lM  Qitk)  +  2AF+(Qi<k,  Qt+hk)]  - 
\  [AF-(Qitk,  Qi+hk)  +  2AF-{Qi+1,k,  Qitk)]  - 


(3.4) 


16 


Where  F  is  the  first-order  numerical  flux  for  the  explicit  Osher's  scheme  given  by 

Fi,k  +  Fi+hk-  fQ'  *  {Fq+-F;}dQ 

J  L/ , 


Fi+i/2,k 


Q.- 
Q, 

± 


whine  I\  =  F+  -f  F~,F*  =  Iqq)  •  and  ■^F±  are  the  corrections  to  obtain 
higher  order  accuracy.  The  Osher  scheme  evaluates  the  flux  assuming  a  shock  tube 
solution  where  Fq  is  piecewise  continuous  and  yields  good  predictions  of  the  flux. 
especially  at  supersonic  Mach  numbers.  For  the  linearization  of  the  left-hand  side  of 
Eq.  3.3,  the  flux  Jacobian  matrices  A,  B  are  evaluated  by  the  Steger-Warming  [Ref. 
28]  flux-vector  splitting. 

Time  accuracy  of  the  implicit  numerical  solution  is  obtained  by  performing 
Newton  iteration  to  convergence  for  each  time  step.  The  approximation  of  Qn+l 
at  each  subiteration  is  the  quantity  Qp.  When  p  >  2,  during  a  given  subiteration, 
Qp  =  Qn+l ,  but  when  p  —  1  and  no  subiteration  is  performed,  then  Qp  =  Qn,  and 
Qp+1  =  Qn+1 .  By  subiterating  to  convergence,  linearization  and  factorization  errors 
are  minimized,  because  the  left-hand  side  of  Eq.  3.3,  can  be  driven  to  zero  at  each 
time  step. 

The  linearization  errors  are  eliminated  by  subiteration  to  convergence.  Typ- 
ically, two  to  three  subiterations  are  sufficient  to  drop  the  residuals  two  orders  of 
magnitude  during  the  Newton  iteration  process. 

High  order  accurate  shock-capturing  schemes  have  some  limitations;  they  may 
select  a  nonphysical  solution,  they  produce  spurious  oscillations  and  they  may  de- 
velop a  nonlinear  instability  in  nonsmooth  and  discontinuous  flow  regions.  More  ap- 
propriate high-order  shock-capturing  schemes,  suitable  for  the  computation  of  weak 
solutions  are  the  TV'D  schemes,  described  in  detail  in  [Ref.  29,  30].  In  the  present 
study,  the  Osher-Chakravarthy  [Ref.  31]  TVD  scheme  is  used.  This  TVD  scheme  has 
flux  limiters  which  impose  constraints  on  the  gradients  of  the  fluxes.  The  flux-limited 
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values  A/     are  computed  from  the  unlimited  fluxes  A/±  as  follows 


&fi+3/2,k  ~  minrnod 

-^■f  t+i/2  k  —  rninmod 
-V^  1  +  1/2,  k  =  minrnod 
m  in  mod 


A/,+_, 


/2.k 


^■Ji+3/2,ki  ^■^■Ji  +  \/2,k 
-^f  1  +  1/2, k'  P&>fi+3/2,k 
^fi+l/2,k>  y^fi-\/2,k 
^fi-l/2,ki  P&fi+l/2,k 


3.5 


where  the  minmod  operator  is  defined  by 


minmod[x,y]  —  sign(x)  x  max[0,min{\x\,ysign{x)}]  (3.6] 

The  viscous  fluxes  Sifk+i/2  are  computed  with  central  differences  as  follows 


■?i,k+l/2  =  S[Qitk+\/2->  {Q<i)i,k+l/2i  Ct,A+l/2J 
Qx,k+\/2  =  n(Qi,k  +  Qi,k-l) 

{Q<:)i,k+i/2  —  Qi,k+i  —  Qi,k 


(3.7; 
(3.s; 

(3.9; 


The  experimental  Reynolds  numbers  based  on  the  chord  length  for  the  test  cases 
examined  are  in  the  range  Rec  %  3.0  x  106  —  5.0  x  106,  and  it  is  expected  that  the  flow 
is  mostly  turbulent.  Transitional  flow  is  expected  to  have  a  small  effect  at  regions  very 
close  to  the  leading  edge.  Present  knowledge  about  boundary  layer  transition  does 
not  enable  computation  and  modeling  of  the  transition  regime.  Therefore,  only  fully 
turbulent  solutions  were  computed.  In  the  present  work,  the  widely  used  two-layer 
Baldwin-Lomax  turbulence  model  was  used.  The  effectiveness  of  other  turbulence 
models,  such  as  the  Johnson-King  model  [Ref.  32]  and  the  RNG  based  algebraic 
model  [Ref.  33]  for  steady  and  unsteady  flows,  was  investigated  in  references  [34]  and 
[35]. 

C.      BOUNDARY  CONDITIONS 

The  solutions  on  the  two  grids  are  computed  separately,  with  the  inner  and 
outer  solutions  communicating  through  the  zonal  interface  boundary.  The  inner  grid 
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surrounds  the  airfoil  and  includes  the  boundary  layer  region  and  the  wake  for  viscous 
solutions.  Inviscid  solutions  are  obtained  by  applying  the  flow  tangency  slip  condition 
where  the  normal  contravariant  velocity  component,  W ',  is  set  equal  to  zero  on  the 
surface.  For  viscous  solutions,  the  nonslip  condition  is  applied  for  the  velocities  on  the 
airfoil  surface.  In  both  cases  the  density  and  pressure  are  obtained  from  the  interior 
grid  points  by  simple  extrapolation.  For  unsteady  solutions,  the  surface  velocity  is 
set  equal  to  the  airfoil  speed  obtained  by  the  prescribed  airfoil  motion  as  follows 

u  =  j((ttz  -  6G),   u>  =  j(Zt(x  -  G6r) 

Unsteady  solutions  for  pitching  and  oscillating  airfoils  are  obtained  by  rotating  the 
inner  grid  only.  Therefore,  only  the  metrics  of  the  inner  grid  must  be  reevaluated  for 
each  time  step. 

At  the  inner  zonal  interface,  the  flow  variables  are  obtained  from  the  interior  of 
the  outer  grid  solution.  Similarly,  the  inner  zonal  boundary  of  the  outer  grid  obtains 
boundary  information  from  the  interior  of  the  inner  grid.  The  inner  and  outer  grid 
radial  lines  are  not  aligned,  in  general.  The  relative  location  of  the  two  grids  with 
respect  to  the  inertial  reference  frame  as  the  inner  grid  rotates  is  computed.  These 
distances  between  neighboring  points  at  the  zonal  interface  are  used  for  a  weighted 
averaging  of  the  conservative  variables. 

All  flows  were  computed  for  subsonic  free-stream  speeds.  At  the  subsonic  inflow 
and  outflow  boundaries  of  the  outer  grid,  the  flow  variables  were  reevaluated  using 
zero-order  Riemann  invariant  extrapolation.  At  the  inflow  boundary,  there  are  three 
incoming  and  one  outgoing  characteristics.  Therefore,  three  variables,  the  density  p, 
the  normal  velocity  w,  and  the  pressure  p,  are  specified  and  the  fourth  variable,  the 
axial  velocity  u  is  extrapolated  from  the  interior.  The  inflow  boundary  conditions  are 
given  by 
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a2    N(l/-K-1) 

Pi     = 


"■»*! 


Poo 


ai  =  l2rl>  (#+_#- 


(3.10 

a,  =  (/?+  +  /?J)/2 

a'!  =  w.x 


Pi  =  l-r1) 


where  /?*.  i?2    are  ^ie  incoming  and  outgoing  Riemann  invariants  given  by 

R\  =  u^  +  20^/(7  -  1).      ftj  =  u2  -  2a2/h  -  1) 

At  the  outflow  boundary  there  are  one  incoming  and  three  outgoing  characteristics. 
I  hereto  re  only  one  quantity,  the  pressure,  is  specified,  while  the  others  are  extrap- 
olated from  the  interior.  For  the  density  and  normal  velocity,  simple  first-order 
extrapolation  is  used,  and  the  axial  outflow  velocity  is  obtained  from  the  zero-ordei 
outgoing  Riemann  invariant.  The  outflow  boundary  conditions  are  given  bv 

P\    =   f>2 

u,  =  flf-2a1/(7-l) 

ai  =  y/lPi/Pi  l:3-'  '  I 

"'l    =   ">2 
Pi    =  P-x 
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were  solved.  For  the  unsteady  flow  solutions,  the  outer  grid  remained  stationary  and 
the  metrics  were  not  reevaluated  at  each  time  step. 
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IV.  RESULTS  AND  DISCUSSION 

The  validity  of  the  zonal  grid  approach  was  first  investigated  for  inviscid  solu- 
tions. An  advantage  of  the  present  approach  is  that  different  grid  densities  may  be 
used  for  the  inner  and  outer  grids.  However,  the  accuracy  and  the  conservative  char- 
acter of  the  solution  for  different  locations  of  the  zonal  interface  and  grid  densities 
must  be  assessed. 

First,  the  accuracy  of  the  computed  results  for  different  inner  and  outer  grid 
densities  was  evaluated.  The  effect  of  the  location  of  the  zonal  interface  relative  to  the 
airfoil  on  the  accuracy  of  the  solution  was  also  investigated.  Then  viscous  solutions  at 
fixed  angles  of  incidence,  up  to  approximately  the  static  stall  angle,  were  computed. 
Finally,  unsteady  flow  responses  to  a  ramp  motion  at  subsonic  free-stream  speed  of 
Moq  =  0.3  and  for  an  oscillation  at  a  free-stream  speed  of  AI^  =  0.6  were  computed. 

A.      STEADY  STATE  SOLUTIONS 
1.      Inviscid  Test  Cases 

Preliminary  test  cases  were  computed  using  coarse  meshes  with  an  inviscid 
solution.  A  two-block  grid  consisting  of  an  81  x  40  point  O-type  inner  grid  and  an 
SI  x  22  point  O-type  outer  grid  was  used  as  a  baseline  grid  for  the  inviscid  solutions. 
Table  4.1  gives  the  inviscid  grids  that  were  tested. 

a.      Case  1.  Baseline  Grid:  81x40  Inner  and  81x22  Outer 

An  inviscid  solution  using  the  baseline  grid  for  subsonic  flow  over  a 
NACA-0012  airfoil  at  M^  —  0.8,  a  =  —.1°  was  obtained.  The  baseline  grid  is  given  in 
figure  4.2.  The  distribution  of  the  computed  surface  pressure  coefficient  is  compared 
with  the  measurements  of  [Ref.  4]  in  Fig.  4.1.  Agreement  with  the  experimental  data 
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TABLE  4.1:   GRID  DENSITIES  OF  THE  INVISCID  SOLUTIONS  COMPUTED 
FOR  A  NACA-0012  AIRFOIL  AT  MM  =  .8  AND  q  =  -0.1. 


Case 

Inner  Grid 

Outer  Grid 

Inner  Grid  Radius 

1 

81  x  40 

81  x  22 

1.5  x  chord 

o 

81  x  40 

41  x  22 

1.5  x  chord 

3 

81  x  40 

81  x  22 

1.0  x  chord 

4 

81  x  40 

81  x  22 

.75  x  chord 

5 

81  x  18 

81  x  19 

OvalGrid 

6 

81  x  20 

41  x  12 

1.5  x  chord 

is  satisfactory  for  an  Euler  solution.  It  can  be  seen  that  with  this  mesh  the  strength 
of  the  shock  is  captured,  but  the  location  is  lagged  by  5  percent  of  the  chord.  The 
baseline  solution  is  conversed  at  2000  iterations. 
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Figure  4.1:  M^  =  0.8  Computed  Solution  With  Baseline  Grid  Compared  to  Exper- 
iment. 
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Figure  4.2:  Inner  Grid  81  x  40  Outer  Grid  is  81  x  22 
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b.      Case  2.   Baseline  Inner  Grid  With  a  41x22  Outer  Grid 

Next  the  effect  of  the  outer  grid  density  is  investigated.  The  grid  for 
Case  2  was  generated  by  starting  with  the  baseline  grid  and  removing  every  other  grid 
point  from  the  outer  grid.  The  resulting  grid  is  given  in  Fig.  4.4.  The  Euler  solution 
is  compared  to  experiment  in  Fig.  4.3.  Overall  the  converged  solution  for  this  case 
agrees  with  the  baseline  solution  except  it  is  noticed  that  the  shock  location  for  the 
upper  and  lower  surface  are  slightly  farther  apart  than  the  baseline  grid  solution. 
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Figure  4.3:   M<x,  —  0.8  Computed  pressure  coefficient  solution  compared  to  experi- 
mental data. 
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Figure  4.4:  Inner  Grid  81  x  40  Outer  Grid  is  41  x  22 
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c.      Case  3.   Overlap  Boundary  Set  to  1  Chord  Length  Away  From 
Body 

In  this  case  the  effect  of  the  overlap  boundary  location  is  studied.  The 
grid  generated  is  given  in  Fig.  4.6.  It  is  seen  that  the  transition  from  the  inner  to 
the  outer  grid  is  not  as  smooth  as  for  the  cases  where  the  boundary  was  located  50 
percent  farther  from  the  airfoil.  At  the  leading  and  trailing  edges,  the  overlapped 
boundary  is  only  a  half  chord  length  away.  In  figure  4.5  the  converged  solution  is 
given.  With  this  grid  very  little  effect  on  the  shock  strength  and  location  is  seen.  The 
pressure  coefficient  before  the  shock  is  slightly  underpredicted  and  after  the  shock  is 
slightly  overpredicted. 
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Figure  4.5:   M^  =  0.8  Computed  Solution  Compared  to  Experiment. 


27 


Figure  4.6:  Inner  Grid  81  x  40  Outer  Grid  is  81  x  22 
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d.      Case  4.     Overlap  Boundary  Set  to   .75  Chord   Length  Away 
From  Airfoil 

In  order  to  further  investigate  the  tendencies  observed  in  Case  3,  an- 
other grid  is  developed  with  the  zonal  interface  closer  to  the  airfoil.  The  overlapped 
region  is  only  a  quarter  chord  away  from  the  leading  and  trailing  edges  of  the  airfoil 
as  seen  in  Fig.  4.8.  The  solution  obtained  is  compared  to  experiment  in  Fig.  4.7. 
The  pressure  coefficient  is  again  slightly  underpredicted  before  the  shock  and  slightly 
overpredicted  after  the  shock.  It  is  also  observed  that  the  zonal  interface  location 
relative  to  the  airfoil  has  little  effect  on  the  solution. 
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Figure  4.7:   M^  =  0.8  Computed  Solution  Compared  to  Experiment. 
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Figure  4.8:  Inner  Grid  81  x  40  Outer  Grid  is  81  x  22 
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e.      Case  5.   Grid  With  Oval  Interface 

The  solution  obtained  for  cases  3  and  4  showed  that  the  location  of 
the  zonal  interface  had  very  little  effect  on  the  computed  solutions.  The  ability  of 
the  zonal  interface  to  pass  flow  discontinuities  was  also  studied.  It  was  known  from 
the  previous  cases  that  the  shock  location  was  near  the  the  mid  chord  point.  The 
inner  grid  for  this  test  case  had  to  be  generated  so  that  the  overlapped  region  was 
very  close  to  the  upper  and  lower  surface  of  the  airfoil.  This  was  accomplished  by 
generating  an  oval  inner  grid  and  an  oval  zonal  interface.  This  grid  is  shown  in  figure 
4.10.  The  computed  surface  pressure  coefficient  distribution  is  shown  in  figure  4.9.  It 
is  in  agreement  with  the  experimental  data  and  with  the  previous  computed  solutions. 
The  computed  flow  quantities,  such  as  density  and  pressure,  showed  that  the  zonal 
approach  used  can  pass  shocks  through  the  zonal  interface.  Figure  4.11  shows  Mach 
contours  which  smoothly  pass  through  the  overlapped  boundary. 
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Figure  4.9:   Moo  =  0.8  Computed  Solution  Compared  to  Experiment. 
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Figure  4.10:   Inner  Grid  81  x  18  Outer  Grid  is  81  x  19 
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Figure  4.11:    Mach  Contours:   Overlapped  Region's  Ability  to  Pass  Flow  Disconti- 
nuities. 


33 


f.      Case  6.   Baseline  Grid  With  Half  the  Radial  Grid  Points  On 
Outer  Grid 

In  the  previous  cases,  the  effect  of  circumferential  resolution  was  stud- 
ied. Next,  the  effect  of  the  radial  resolution  on  the  computed  solutions  is  investigated. 
Figure  4.13  shows  the  grid  generated  to  test  these  effects.  The  computed  surface 
pressure  coefficient  distribution  (in  figure  4.12)  was  comparable  with  the  solutions 
obtained  with  denser  outer  grids.  Computations  with  an  even  coarser  grid,  e.g.,  a 
41  x  11  point  grid,  predicted  the  shock  location  even  further  downstream  due  to  lack 
of  stream  wise  resolution. 
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Figure  4.12:   A/^  =  0.8  Computed  Solution  Compared  to  Experiment. 
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Figure  4.13:  Inner  Grid  81  x  40  Outer  Grid  is  41  x  12 
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2.      Viscous  Test  Case 

The  Euler  solutions  presented  in  the  previous  section  were  used  to  study 
the  effects  of  different  grid  densities  on  the  computed  solutions.  Having  confidence 
in  the  predictions  of  the  code,  the  next  step  was  to  generate  a  viscous  grid  with  the 
airfoil's  quarter  chord  point  at  the  center  of  the  inner  grid. 

Viscous,  subsonic  flow  solutions  were  obtained  at  the  following  fixed  angles 
of  attack:  3.27°,  4.97°,  6.69°,  8.38°.  9.27°,  10. 12°,  10. 99°  and  11.90°.  The  flow  condi- 
tions of  the  measurements  reported  in  [Ref.  21]  were  used,  e.g.,  M^  =  0.3,  Re  = 
4.0  x  106.  The  spacing  of  the  grid  was  set  to  0.0005  for  the  first  grid  point  above 
the  surface.  The  quarter  chord  point  of  the  airfoil  was  set  at  the  center  of  rotation 
so  the  airfoil  could  be  ramped  about  the  quarter  chord  point.  These  solutions  were 
obtained  on  a  1S1  x  56  point  viscous  inner  grid  and  a  181  x  26  point  inviscid  outer 
grid  shown  in  figure  4.14. 

Solutions  were  also  computed  on  a  grid  with  half  the  streamwise  resolution. 
e.g..  a  91  x56  point  grid.  The  computed  surface  pressure  coefficient  distributions  using 
the  181  x  56  inner  grid  and  181  x  26  outer  grid,  are  compared  to  experimental  results 
in  figures  4.15  through  4.23. 

Solutions  for  fixed  angles  of  incidence  were  obtained  by  two  methods.  First 
by  rotating  the  inner  grid  to  the  specified  angle  of  incidence  and  setting  the  oncoming 
flow  to  zero  degrees.  Second  by  rotating  the  flow  to  the  angle  of  incidence  and 
leaving  the  inner  grid  at  zero  angle  relative  to  the  outer  grid.  The  computed  pressure 
coefficients  and  boundary  layers  were  the  same  for  both  cases. 

For  the  low  Mach  number  viscous  solutions,  no  flux  limiting  was  applied. 
It  is  seen  that  the  computed  results  closely  agree  with  the  experimental  data.  At  the 
higher  angles  of  incidence  the  suction  peak  is  not  exactly  captured.  This  is  probably 
due  to  lack  of  grid  resolution  at  the  leading  edge  of  the  airfoil. 
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Figure  4.14:    A  NACA-0012  Airfoil  Grid  Centered  at  the  Quarter  Chord  Point, 
181x58  Inner  Grid  and  181x26  Outer. 
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Figure  4.15:   Viscous  Computed  Solution  Compared  to  Experiment  for  a  =  3.27°. 
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Figure  4.16:   Viscous  Computed  Solution  Compared  to  Experiment  for  a  —  4.97°. 
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Figure  4.17:   Viscous  Computed  Solution  Compared  to  Experiment  for  a  —  6.69°. 
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Figure  4.18:  Viscous  Computed  Solution  Compared  to  Experiment  for  a  =  7.54c 
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Figure  4.19:   Viscous  Computed  Solution  Compared  to  Experiment  for  a  =  S.38°. 
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Figure  4.20:  Viscous  Computed  Solution  Compared  to  Experiment  for  a  —  9.27°. 
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Figure  4.21:   Viscous  Computed  Solution  Compared  to  Experiment  for  a  =  10.12°. 


o 


c 

o 

it 

0) 

o 
O 

</> 

0) 


a=10.99  deg.,  M=0.30,  Re=2700000 

-i 1 1 1 1 r 


-0.2        0 


0.2       0.4        0.6       0.8 
Axial  Location  x/c 


Figure  4.22:   Viscous  Computed  Solution  Compared  to  Experiment  for  a  =  10.99". 
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Figure  4.23:   Viscous  Computed  Solution  Compared  to  Experiment  for  a  =  11.99°. 
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B.      RAMP  MOTION  SOLUTION 

The  unsteady  solution  for  a  ramp  motion  from  a  —  0°  to  a  =  30°  at  M,^  = 
0.3,  Re  =  2.7  x  106  and  reduced  frequency  ,  k  =  0.0127  was  obtained  on  both  a 
91  x  56  point  inner  grid  and  a  181  x  56  point  inner  grid.  The  pitch  rate  for  the  ramp 
motion  k  is  defined  as  k  =  ctc/'lU^.  The  computed  lift  response  is  compared  with 
the  experimental  measurements  of  [Ref.  6]  in  Figure  4.24. 
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Figure  4.24:     Comparison  of  the  Measured  and  Computed  Lift  for  the  Ramp  Motion 

Both  the  coarse  and  fine  grid  solutions  closely  predict  the  measured  lift.  How- 
ever, at  the  higher  angles  of  attack,  the  finer  grid  gives  higher  lift.  The  computed 
surface  pressure  coefficient  distributions  at  several  angles  of  incidence  are  compared  in 
Figures  4.25  -  4.38.  Experimental  surface  pressure  coefficients  were  available  for  the 
following  angles  of  incidence:  2.94°, 5.84°, 8.91°,  11.76°,  15.5°,  and  they  are  displayed 
as  diamonds  in  the  figures.  The  computed  surface  pressure  coefficient  distribution 
is  in  good  agreement  with  the  measured  data  over  the  entire  incidence  range.  The 
computed  flowfield  at  the  maximum  angle  of  incidence,  q  =  15.5°,  is  mostly  attached. 
A  small  separated  flow  region  exists  at  the  trailing  e  .ge  region  only.  At  a  higher  angle 
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of  attack,  a  %  17.0°.  the  computed  solution  shows  the  development  of  the  dynamic 
stall  vortex  in  the  leading  edge  region. 
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Figure  4.25:     Computed  ramp  solution  at  a  =  1.86°. 
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Figure  4.26:     Computed  ramp  solution  compared  to  experimental  data  at  a  =  2.94c 
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Figure  4.27:     Computed  ramp  solution  at  a  =  4.14°. 
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Figure  4.28:     Computed  ramp  solution  at  a  =  4.87' 
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Figure  4.29:     Computed  ramp  solution  compared  to  experimental  data  at  q  =  5.85°. 
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Figure  4.30:     Computed  ramp  solution  at  a  —  6.72°. 
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Figure  4.31:     Computed  ramp  solution  at  a  =  8.02°. 
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Figure  4.32:     Computed  ramp  solution  compared  to  experimental  data  at  a  =  S.91w. 
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Figure  4.33:     Computed  ramp  solution  at  a  =  9.85°. 
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Figure  4.34:     Computed  ramp  solution  at  a  =  10.80°. 
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Figure    4.35: 

a  =  11.77°. 
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Figure  4.36:     Computed  ramp  solution  at  a  =  12.84°. 
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Figure  4.37:     Computed  ramp  solution  at  a  =  13.S90. 
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Figure    4.38: 

q  =  15.55°. 


Computed    ramp   solution   compared   to   experimental   data   at 
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1.      Boundary  Layer  Comparisons  For  Ramp  Motion 

Early  on  in  the  development  of  the  software,  the  computed  pressure  coef- 
ficients agreed  quite  well  with  the  experimental  data.  As  test  cases  at  higher  angles- 
of-attack  were  tested,  it  was  discovered  that  the  flow  was  not  separating  from  the 
trailing  edge  of  the  airfoil.  The  problem  was  discovered  by  investigating  the  boundary 
layer  profiles.  A  problem  in  the  turbulence  model  was  discovered  and  easily  fixed. 
The  lesson  here  is  that  although  the  computed  and  experimental  pressure  coefficients 
agree  quite  well  it  is  important  to  check  that  the  boundary  layer  profiles  are  reason- 
able. In  figures  4.39  through  4.48  the  computed  boundary  layers  are  compared  to 
boundary  layers  computed  by  the  interactive  viscous  inviscid  boundary  layer  method 
of  reference  [22]. 

The  comparisons  are  made  at  x/c  =  .5  and  at  x/c  —  .9  for  angles-of- 
attack  ranging  from  2.94°  to  15.5°.  The  computed  boundary  layer  profiles  compare 
quite  well  up  to  15.5°.  At  x/c  =  .9  for  the  angle-of-attack  of  15.5°,  the  comparisons 
diverge.  This  happens  because  of  the  different  turbulence  models  used.  At  this  angie- 
of  attack,  the  flow  is  separating  at  the  trailing  edge  so  the  turbulence  models  used 
become  very  important. 
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Figure  4.39:  Comparison  of  computed  boundary  layer  with  an  interactive  boundary 
layer  program  at  the  90%  chord  for  a  =  2.94  degrees. 
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Figure  4.40:  Comparison  of  computed  boundary  layer  with  an  interactive  boundary 
layer  program  at  the  90%  chord  for  a  =  2.94  degrees. 
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Figure  4.41:  Comparison  of  computed  boundary  layer  with  an  interactive  boundary 
layer  program  at  the  50%  chord  for  a  —  5.85  degrees. 
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Figure  4.42:  Comparison  of  computed  boundary  layer  with  an  interactive  boundary 
layer  program  at  the  90%  chord  for  a  =  5.85  degrees. 
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Figure  4.43:  Comparison  of  computed  boundary  layer  with  an  interactive  boundary 
layer  program  at  the  50%  chord  for  a  =  8.91  degrees. 
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Figure  4.44:  Comparison  of  computed  boundary  layer  with  an  interactive  boundary 
layer  program  at  the  90%  chord  for  a  =  8.91  degrees. 
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Figure  4.45:  Comparison  of  computed  boundary  layer  with  an  interactive  boundary 
layer  program  at  the  50%  chord  for  a  =  11.77  degrees. 
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Figure  4.46:  Comparison  of  computed  boundary  layer  with  an  interactive  boundary 
layer  program  at  the  90%  chord  for  a  =  11.77  degrees. 
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Figure  4.47:  Comparison  of  computed  boundary  layer  with  an  interactive  boundary 
layer  program  at  the  50%  chord  for  a  =  15.55  degrees. 
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Figure  4.48:  Comparison  of  computed  boundary  layer  with  an  interactive  boundary 
layer  program  at  the  90%  chord  for  a  =  15.55  degrees. 
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2.      Computed  Ramp  Flow  Details 

In  this  section,  flow  features  of  the  computed  ramp  motion  arc  shown. 

ir-«'  features  are  demonstrated  by  the  density  contour  lines.  Mach  contour  lines. 

rticity  contour  lines  and  mass-flux  contours.  The  contour  lines  represent  areas  ol 
i  he  flowfield,  for  which  the  parametei  of  interest,  remains  constant.  Clustered  contoui 
lines  represent  areas  where  the  parameter  is  rapidly  changing. 

In  figures.  4.49  through  4.64.  density,  Mach  number,  vorticity  and  mass 
II u.\  contours  are  displayed  for  a  —  1.86"  to  o  =  17.50".   In  figures  4.65  through  4.1  10 
the  Mach  number  contours  are  replaced  by  pressure  contours.  This  was  done  because 
i  he  strong  gradients,  caused  by  the  vortices,  tend  to  not  give  clear  Mach  numbei 
information. 

In  figure  4.63.  a  =  16.50°.  the  beginning  of  a  vortex  is  visible  near  the 
leading  edge.  At  a  =  17.00°.  in  figure  4.64  the  vortex  has  moved  5/1  of  chord 
downstream,  and  another  vortex  is  starting  at  the  leading  edge.  At  n  =  10.00 
i  he  first  two  vortices  merge  into  one  stronger  vortex  located  near  the  509c  chord 
point .  Two  smaller  vortices  are  also  clearly  visible  in  the  mass-flux  contours  of  figure 
1. 68.  The  original  vortex  now  detaches  from  the  surface  of  the  airfoil  at  19.50  .  At 
(i  =  20.50°,  a  second  vortex  is  shed  from  the  leading  edge  area,  and  a  the  vortex  shown 
.it  10. 5°  has  moved  to  the  trailing  edge.  At  a  =  22.50°,  the  trailing  edge  vortex  has 
been  sited  downstream.  The  whole  process  seems  to  repeat  itself  at  a  faster  rate,  as 
can  be  seen  in  the  remaining  figures  as  the  airfoil  ramps  up  to  30".  The  dynamic  stall 
phenomenon  is  observed  in  terms  of  density,  pressure  and  vorticity.  The  nose-down 
pitching  moment  is  caused  by  the  vortex  sitting  at  the  trailing  edge.  It  is  important 
l  hi'  the  vortices  do  not  dissipate  and  do  not  get  distorted  as  they  pass  through  the 
zonal  interface. 


01 


Next  the  solution  was  continued  to  angles  beyond  stall.  Traces  of  alternate 
vortex  shedding  from  the  trailing  edge  can  be  seen  at  q  %  27c'.    At  high  angles  ol 


incidence  the  alternating  vortex  shedding  is  very  well  demonstrated  in  figures  I.1) 
llirough  4.110.  These  computed  solutions  are  in  general  agreement  with  the  finding 
of  the  experimental  investigations  of  Chandrasekahara  et  al.  Also  it  can  be  seen  tha 
i  he  zonal  method  developed  is  capable  of  computing  unsteady  flows  at  very  higl 
,111'jjes-of-attack  showing  at  least  qualitative  agreement  with  the  experiment. 
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igure  4.49:    Ramp  Motion  Flow  Details,  1 
=  1.86°. 


Moo  =  .3,  Jfc  =  .0127,  #e  =  2.7  x  106, 
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Figure  4.50:    Ramp  Motion  Flow  Details,  A/c 
a  =  2.94°. 


=  .3,  k  =  .0127,  Re  =  2.7  x  106. 
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Figure  4.51:    Ramp  Motion  Flow  Details,  M0 

a  =  4.14°. 


=  .3,  k  =  .0127,  Re  =  2.7  x  10* 
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Figure  4.52:    Ramp  Motion  Flow  Details,  MM  =  .3,  fc  =  .0127,  #e  =  2.7  x  106, 
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Figure  4.53:    Ramp  Motion  Flow  Details,  M0 
a  =  6.72°. 


=  .3,  k  =  .0127,  7?e  =  2.7  x  10f 
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Figure  4.54:    Ramp  Motion  Flow  Details,  M0 

a  =  8.02°. 


=  .3,  k  =  .0127,  Re  =  2.7  x  10* 
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Figure  4.55:    Ramp  Motion  Flow  Details,  M0 
a  =  8.91°. 


=  .3,  k  =  .0127,  Re  =  2.7  x  106. 
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Figure  4.56:    Ramp  Motion  Flow  Details,  M0 
a  =  9.85°. 


=  .3,  k  =  .0127,  Re  =  2.7  x  106, 
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Figure  4.57:    Ramp  Motion  Flow  Details,  A/c 
a  =  10.80°. 
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Figure  4.58:    Ramp  Motion  Flow  Details,  M^  =  .3,  k  =  .0127,  i?e  =  2.7  x  106. 
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Figure  4.59:    Ramp  Motion  Flow  Details,  A/0 
a  =  12.84°. 
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Figure  4.60: 

a  =  13.89°. 


Ramp  Motion  Flow  Details,  A/^  =  .3,  k  =  .0127,  fie  =  2.7  x  106. 
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Figure  4.61:    Ramp  Motion  Flow  Details,  M0 
a  =  15.55°. 
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Figure  4.62: 

a  =  16.00°. 


Ramp  Motion  Flow  Details,  M^  =  .3,  k  =  .0127,  Re  =  2.7  x  106. 
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Figure  4.63:    Ramp  Motion  Flow  Details,  Mc 
a  =  16.50°. 


=  .3,  k  =  .0127,  Re  =  2.7  x  106, 
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Figure  4.64:    Ramp  Motion  Flow  Details,  Mc 
a  =  17.00°. 


=  .3,  k  =  .0127,  fle  =  2.7  x  106, 
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Figure  4.65: 

a  =  17.50°. 


Ramp  Motion  Flow  Details,  M0 


.3,  k  =  .0127,  #e  =  2.7  x  106. 
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Figure  4.66:    Ramp  Motion  Flow  Details,  Mc 
a  =  18.00°. 


=  .3,  k  =  .0127,  fle  =  2.7  x  10( 
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Figure  4.67: 

q  =  18.50°. 


Ramp  Motion  Flow  Details,  M^  =  .3,  k  =  .0127,  Re  =  2.7  x  106. 
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Figure  4.68:    Ramp  Motion  Flow  Details,  M0 
q  =  19.00°. 


=  .3,  k  =  .0127,  Re  =  2.7  x  106, 
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Figure  4.69:    Ramp  Motion  Flow  Details,  Mc 
a  =  19.50°. 


=  .3,  k  =  .0127,  Re  =  2.7  x  106. 
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Figure  4.70:    Ramp  Motion  Flow  Details,  Mc 
q  =  20.00°. 


=  .3,  A:  =  .0127,  Re  =  2.7  x  106. 
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Figure  4.71:    Ramp  Motion  Flow  Details,  M0 
a  =  20.50°. 


=  .3,  k  =  .0127,  #e  =  2.7  x  106, 
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Figure  4.72:    Ramp  Motion  Flow  Details,  A/0 
a  =  21.00°. 


=  .3,  k  =  .0127,  Re  =  2.7  x  106, 
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Figure  4.73:    Ramp  Motion  Flow  Details,  M0 
a  =  21.50°. 


=  .3,  k  =  .0127,  #e  =  2.7  x  106. 
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Figure  4.74: 

a  =  22.00°. 


Ramp  Motion  Flow  Details,  M^  =  .3,  k  =  .0127,  Re  =  2.7  x  106, 
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Figure  4.75:    Ramp  Motion  Flow  Details,  M0 
a  =  22.50°. 


=  .3,  k  =  .0127,  #e  =  2.7  x  106, 
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Figure  4.76:    Ramp  Motion  Flow  Details,  M0 
a  =  23.00°. 


-  .3,  k  =  .0127,  7?e  =  2.7  x  106, 
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Figure  4.77:    Ramp  Motion  Flow  Details,  Mc 
q  =  23.50°. 


=  .3,  fc  =  .0127,  Re  =  2.7  x  106, 
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Figure  4.78: 
a  =  24.00°. 


Ramp  Motion  Flow  Details,  MM  =  .3,  fc  =  .0127,  Re  =  2.7  x  106, 
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Figure  4.79:    Ramp  Motion  Flow  Details,  Mc 
a  =  24.50°. 


=  .3,  k  =  .0127,  Re  =  2.7  x  106. 
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Figure  4.80:    Ramp  Motion  Flow  Details,  M0 
a  =  25.00°. 


=  .3,  k  =  .0127,  Re  =  2.7  x  106, 
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Figure  4.81:    Ramp  Motion  Flow  Details,  M0 
a  =  25.50°. 


=  .3,  k  =  .0127,  Re  =  2.7  x  106. 
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Figure  4.82:    Ramp  Motion  Flow  Details,  MM  =  .3,  A;  =  .0127,  Re 
a  =  26.00°. 
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Figure  4.83:    Ramp  Motion  Flow  Details.  Mz 

<i  =  26.50°. 
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Figure  4.84:    Ramp  Motion  Flow  Details,  A/0 

n  =  27.00". 


=  .3,  k  =  .0127,  #e  =  2.7  x  10* 
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Figure  4.85:    Ramp  Motion  Flow  Details,  M0 
a  =  27.50°. 


=  .3,  k  =  .0127,  Re  =  2.7  x  10fc 
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Figure  4.86:    Ramp  Motion  Flow  Details,  Mc 
a  =  28.00°. 


=  .3,  k  =  .0127,  #e  =  2.7  x  10* 
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Figure  4.87:    Ramp  Motion  Flow  Details.  Mc 
a  =  28.50°. 
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Figure  4.88:    Ramp  Motion  Flow  Details.  M0 
m  =  29.00°. 
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Figure  4.89: 
n  =  29.50°. 


Ramp  Motion  Flow  Details.  M*,  =  .3,  k  =  .0127,  Re  =  2.7  x  10* 
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Figure  4.90:    Ramp  Motion  Flow  Details.  Mc 
n  =  :?0.00°. 
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Figure  4.91:    Ramp  Motion  Flow  Details,  M0 
<,  =  31.00°. 


=  .3,  k  =  .0127,  /?e  =  2.7  x  10b. 
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Figure  4.92:    Ramp  Motion  Flow  Details.  .V/c 
n  =  32.00°. 


=  .3,  k  =  .0127,  #e  =  2.7  x  106. 
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Figure  4.93:    Ramp  Motion  Flow  Details.  Mc 
n  =  33.00°. 


=  .3,  k  =  .0127,  i?e  =  2.7  x  10b. 
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Figure  4.94:    Ramp  Motion  Flow  Details,  MQ 
a  =  34.00°. 


=  .3,  k  =  .0127,  #e  -  2.7  x  106, 
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Figure  4.95:    Ramp  Motion  Flow  Details,  M0 
a  =  35.00°. 


=  .3,  k  =  .0127,  Re  =  2.7  x  10s 
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Figure  4.96: 

q  =  36.00°. 


Ramp  Motion  Flow  Details,  M^  =  .3,  k  =  .0127,  Re  =  2.7  x  106, 
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Figure  4.97:    Ramp  Motion  Flow  Details.  Mc 
(i  =  37.00°. 


=  .3,  fc  =  .0127,  i?e  =  2.7  x  106. 
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Figure  4.98: 
a  =  38.00°. 


Ramp  Motion  Flow  Details,  M^  =  .3,  k  =  .0127,  i?e  =  2.7  x  106, 
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Figure  4.99:    Ramp  Motion  Flow  Details.  A/0 
n  =  :W.00°. 


=  .3,  fc  =  .0127.  #e  =  2.7  x  lO11 
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Figure  4.100: 

m  =  10.00". 


Ramp  Motion  Flow  Details.  M^  =  .3,  k  =  .0127,  Re  =  2.7  x  10*' 
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Figure  4.101:    Ramp  Motion  Flow  Details,  M0 
a  =  41.00°. 


=  .3,  k  =  .0127,  Re  =  2.7  x  10f 
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Figure  4.102:    Ramp  Motion  Flow  Details.  M^  =  .3,  k  =  .0127,  Re  =  2.7  x  10' 

M     =      1:2.00°. 
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Figure  4.103:    Ramp  Motion  Flow  Details.  M0 

n    =    13.00". 
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Figure  4.104:    Ramp  Motion  Flow  Details,  M0 

n  =  14.00°. 
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Figure  4.105:    Ramp  Motion  Flow  Details.  M0 

M    =     to. 00". 
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Figure  4.106:    Ramp  Motion  Flow  Details,  M^  =  .3,  k  =  .0127,  Re  =  2.7  x  10e 

q  =  46.00°. 
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Figure  4.107:    Ramp  Motion  Flow  Detail 
..  =  t7.00L\ 


.U^  =  .3.  k  =  .0127,  Re  =  2.7  x  10b 
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Figure  4.108:    Ramp  Motion  Flow  Details.  M«,  -  .3,  k  =  .0127,  Re  =  2.7  x  10fc 
n  =  tS.OO0. 


118 


on 

3 
(/) 
00 

UJ 


>- 

o 

I— 
an 
O 
> 


Figure  4.109:    Ramp  Motion  Flow  Details.  Ma 
n  =  -19.00°. 


=  .3,  A;  =  .0127,  i?e  =  2.7  x  106 
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Figure  4.110:    Ramp  Motion  Flow  Details,  Mc 
m  =  -)0.00°. 


=  .3,  k  =  .0127,  Re  =  2.7  x  10e 
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C  .      OSCILLATORY  MOTION  SOLUTION 

1  he  unsteady  solution  for  a  periodic  oscillatory  motion,  given  by  a{t)  =  4.86  + 
_'.  I  1  sin[vt),  at  Moo  =  0.6.  Rec  =  4.S  x  10".  with  a  reduced  frequency  of  k  =  O.K. 
A-as  also  obtained.  Here  the  reduced  frequency  is  defined  as  k  =  uje/Urx,.  The  flow  foi 
■  his  motion  is  initially  purely  subsonic;  but.  as  the  angle  of  attack  increases  to  about 
nf  / )  ~  o° .  supersonic  flow  conditions  are  encountered  at  the  leading  edge  region  and  i\ 
i  ransonic  shock  forms.  This  shock  is  present  during  the  upstroke  until  the  maximum 
angle  of  attack  is  reached  and  during  the  downstroke  up  to  about  a(t)  %  5.0°.  Th<' 
computed  and  measured  lift  and  pitching  moment  response  are  compared  in  Figs. 
I.I  I  1  and  4.112.  respectively. 
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Figure  4.111:    Comparison  of  Measured  and  Computed  Lift  for  Oscillatory  Tesi 
Moi  ion 

The  computed  lift  and  pitching  moment  coefficients  are  in  close  agreement  with 
the  measured  values.  The  computed  surface  pressure  distribution  is  compared  with 
the  measurements  of  reference  [6]  for  two  angles  during  the  upstroke  and  two  angle- 
during  the  downstroke  in  Figs.  4.113  through  4.116  The  computed  surface  pressure  is 
in  better  agreement  with  the  measurements  at  the  lower  angles  of  incidence  (o  =  5.95 ' 
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up  and  q  =  5.11°  clown).  At  higher  incidences  (Fig.  4.114  and  4.115),  the  agreement 
lUMei'iorates  in  the  region  around  the  shock.  The  global  view  of  the  computed  density 
Held  shows  that  the  density  contours  smoothly  cross  the  zonal  interface  for  the  case 
\\  here  a  shock  exists. 

a(t)=4.86+2.44sin(wt),  M=0.60,  k=0.16,  Re=4800000 
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Figure   4.112:     Comparison  of  Measured  and  Computed  Moment  Coefficient  foi 
Oscillatory  Test  Case 
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Figure  4.113:  Comparison  of  the  Measured  and  Computed  Unsteady  Surface  Pres- 
sure Coefficient  of  Oscillatory  Motion,  o  =  5.95°  upstroke. 
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Figure  4.114:   Comparison  of  the  Measured  and  Computed  Unsteady  Surface  Pro- 
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Coefficient  of  Oscillatory  Motion,  a  =  6.97°  upstroke. 
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Figure  4.115:   Comparison  of  the  Measured  and  Computed  Unsteady  Surface  Pres- 
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Coefficient  of  Oscillatory  Motion,  o  =  6.57°  downstroke. 


a=5.1l  deg.  down,  M=0.60,  k=0.16,  Re=4800000 
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Figure  4.116:    Comparison  of  the  Measured  and  Computed  Unsteady  Surface  Pres- 
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e  Coefficient  of  Oscillatory  Motion, 


.95°  downstroke. 
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V.  CONCLUSIONS 

A  solution  procedure  suitable  for  steady  and  unsteady  compressible  flow  solu- 
tions using  zonal  overlapped  grids  was  developed.  Simple  weighted  averaging  \va> 
u-ed  at  the  overlapped  zonal  interfaces.  Steady  and  unsteady,  inviscid  and  viscoii- 
lluw  solutions  for  subsonic  and  transonic  flows  over  airfoils  were  presented  to  validate 
i  he  zonal  grid  approach. 

The  inviscid  solutions  presented  show  the  overlapped  zonal  interface's  ability  to 
pass  flow  properties  without  distortion.  It  was  found  that  for  the  inviscid  test  cases 
I  he  location  of  the  zonal  interface  is  not  important.  In  fact,  as  the  zonal  interface 
moved  closer  to  the  airfoil,  while  keeping  the  number  of  grid  points  constant,  the 
pressure  coefficient  prediction  actually  improved.  This  is  due  to  more  grid  point- 
being  clustered  close  to  the  airfoil.  This  test  case  had  strong  shocks  on  the  upper  and 
lower  surface  of  the  airfoil,  and  as  the  grid  points  were  moved  closer  to  the  airfoil, 
i  lie  computed  solution  tended  to  have  bigger  oscillations  near  the  shock. 

1  he  steady  viscous  test  cases  showed  that  using  the  Baldwin-Lomax  turbulent  <■ 
model  with  the  present  approach  gave  accurate  results  for  attached  and  mildly  sep- 
arated flow  over  stationary  airfoils.  This  case  demonstrates  one  of  the  advantage- 
of  the  present  approach,  specifically,  that  solving  the  inviscid  equations  on  the  outei 
Lirid  and  the  viscous  equations  on  the  inner  grid  gives  good  results.  The  flow  variables 
<  re  also  passed  smoothly  through  the  zonal  interface. 

The  ramp  case  again  showed  good  agreement  with  experimental  data.  With 
this  rase  another  advantage  of  the  present  approach  was  displayed.  The  inner  grid 
was  rotated  with  the  airfoil  to  the  new  angles  ot  attack,  i  he  high  order  accural' 
-(heme  enabled  the  software  to  convect  vortices,  and  at  high  angles,  these  vortices 
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convected  up  to  3  chord  lengths. 

The  final  test  case  was  an  oscillatory  one.  For  this  test  case,  supersonic  flow  was 
encountered  on  the  upstroke  at  the  leading  edge  region  which  produced  a  transonic 
shock.  The  computed  lift  and  pitching  moment  coefficients  are  in  close  agreement 
with  the  measured  values.  At  the  higher  angles  of  attack  the  agreement  with  measured 
data  deteriorated  in  regions  near  the  shock. 

The  following  are  some  recommendations  based  on  this  study. 

1.  It  has  been  shown  that  the  zonal  grid  approach  can  be  used  to  study  unsteady 
viscous  flows.  A  systematic  comparison  with  other  codes  should  be  conducted 
in  order  to  quantify  the  efficiency  of  the  present  approach. 

2.  Presently,  there  exists  no  easy  way  to  generate  zonal  grids.  In  this  study  the 
grids  were  generated  separately  and  then  refined  through  several  other  pro- 
grams. The  development  of  a  new  or  the  modification  of  an  existing  software 
package  such  as  GRAPE  is  recommended.  The  user  should  be  able  to  specify, 
along  with  all  the  usual  information,  airfoil  shape,  the  location  of  the  overlapped 
region,  the  number  of  overlapped  cells  and  the  size  of  the  outer  grid. 

3.  An  advantage  of  the  present  approach  is  that  it  can  be  extended  to  multiple 
inner  and  outer  grids.  A  flow  solver  needs  to  developed  to  take  advantage  of 
this  so  that  effects  of  oscillating  bodies  in  relative  motion  can  be  studied. 
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